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It was recently shown that monotone gene action, i.e., order-preservation between 
allele content and corresponding genotypic values in the mapping from genotypes to 
phenotypes, is a prerequisite for achieving a predictable parent-offspring relationship 
across the whole allele frequency spectrum. Here we test the consequential prediction 
that the design principles underlying gene regulatory networks are likely to generate 
highly monotone genotype-phenotype maps. To this end we present two measures 
of the monotonicity of a genotype-phenotype map, one based on allele substitution 
effects, and the other based on isotonic regression. We apply these measures to 
genotype-phenotype maps emerging from simulations of 1881 different 3-gene regulatory 
networks. We confirm that in general, genotype-phenotype maps are indeed highly 
monotonic across network types. However, regulatory motifs involving incoherent 
feedforward or positive feedback, as well as pleiotropy in the mapping between genotypes 
and gene regulatory parameters, are clearly predisposed for generating non-monotonicity. 
We present analytical results confirming these deep connections between molecular 
regulatory architecture and monotonicity properties of the genotype-phenotype map. 
These connections seem to be beyond reach by the classical distinction between additive 
and non-additive gene action. 
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INTRODUCTION 

Quantitative genetics is the major theoretical foundation for 
genetic studies in production biology, evolutionary biology, and 
biomedicine. A core concept in quantitative genetics is the geno- 
typic value, the mean observed phenotype for a given genotype. 
It constitutes the basis for the genotype-to-phenotype (GP) map 
concept. The shape of a given GP map is typically described 
by the classical gene action terms: additivity, dominance, and 
epistasis. Together with genotype frequencies in a given pop- 
ulation, the GP map is the basis for decomposing observed 
phenotypic variance into environmental variance and genetic 
variance components including additive variance, dominance 
variance and epistatic variance. This provides the basis for a 
very successful theory when it comes to predicting selection 
response and breeding values (Falconer and Mackay, 1996; Lynch 
and Walsh, 1998) and more recent statistical genetics meth- 
ods for mapping Quantitative Trait Loci (QTL) (Neale et al., 
2008). Quantitative genetics thus provides a mature machin- 
ery for predicting the population level consequences of a given 
GP map, but in order to understand several generic genetic 
phenomena there is a stated need for new tools for disclos- 
ing how the shape of the GP map is determined by underly- 
ing biology (Jaeger et al., 2012; Moore, 2012; Gjuvsland et al., 
2013). 

One such phenomenon is the resemblance between parents 
and offspring. An explanation in quantitative genetic terms is 
that the additive variance (V^) makes up a substantial part of 



the phenotypic (Vp) and genetic variance {Vq). Hill et al. (2008) 
showed that in populations with extreme allele frequencies, high 
Va/Vg ratios will arise regardless of the shape of the GP map. 
However, for populations with intermediate allele frequencies a 
much wider range of Va/Vq ratios is observed (Wang et al., 2013). 
In such populations, high VaIVg ratios cannot be fully accounted 
for without considering properties of the GP map. Gjuvsland 
et al. (2011) showed that a key feature of GP maps that give 
high ratios of additive to genotypic variance (Va/Vg), * s a mono- 
tone (or order-preserving) relation between gene content (the 
number of alleles of a given type) and phenotype. This led to 
the hypothesis that the regulatory circuitry of sexually reproduc- 
ing organisms predominantly predisposes for highly monotone 
genotype-phenotype maps. 

Here we address the above hypothesis by a two-step approach. 
First we provide methods and software tools for measuring 
monotonicity of generic GP maps (i.e., sets of genotypic values). 
Then we use these tools on the data generated by an extensive 
simulation study of a broad collection of gene regulatory net- 
work models. In these network models the steady state expression 
levels serve as phenotypes and genetic variation is introduced 
through parameters describing maximal production rates and 
the shape of the gene regulation function. Such causally cohe- 
sive genotype-phenotype (cGP) models [see Gjuvsland et al. (2013) 
and references therein] allow us to identify relationships between 
regulatory network architecture and properties of the resulting 
GP maps. 
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Our results confirm the prediction that the GP maps arising 
from a wide range of gene regulatory network motifs are in gen- 
eral highly monotone. In addition we show through numerical 
as well as mathematical analysis that regulatory motifs involv- 
ing incoherent feed-forward or positive feedback stand out in 
their capacity to generate non-monotonicity. These relationships 
between molecular regulatory architecture and properties of the 
genotype-phenotype map — of substantial relevance to functional 
genomics in general — are beyond reach by the standard distinc- 
tion between additive and non-additive gene action. 

Our approach can be applied to cGP models of a wide range of 
biological systems at any level of model complexity. It opens for 
a systematic study of the monotonicity properties of molecular 
regulatory structures underlying the whole spectrum of physio- 
logical regulation. This suggests that the concept of monotonicity 
of GP maps can be used to build theory about heredity phrased 
in terms of molecular mechanism, something which standard 
genetic concepts and approaches appear to be incapable of. 

MODELS AND METHODS 

BACKGROUND ON MONOTONICITY OF GP MAPS 

To ease understanding we provide a brief recapitulation of 
the concept of monotonicity (or order-preservation) in GP 
maps introduced in (Gjuvsland et al, 2011). We consider a 
diploid genetic model with N biallelic loci (alleles indexed 1 
and 2) underlying a quantitative phenotype. A genotype at a 
single locus k is denoted by g k € {11, 12,22}. In the case of 
two loci k and 1 there are 9 possible genotypes gti = gkgi € 
{1111, 1112, 1122, 1211, 2212, 2222}. The general N loci 
genotype space F contains 3 N genotypes gig2 ■ ■ ■ gN (in con- 
densed notation gi : jv) constructed by concatenating single locus 
genotypes, V = [g lg2 ■■■g N \gk& {11, 12, 22}, k= 1,2, ... , N}. 



For any locus k, the genotypic background, i.e., the allele com- 
position of all loci except k, is g (k) = gig 2 ■ ■ -gk-igk+i ■ ■ -gN = 
gi-.k-\gk+v.N- For example, if N = 4 then g (2 ' = 112212 
means that the genotypes of locus 1, 3, and 4 are 11, 
22, and 12, respectively. We use the straightforward nota- 
tion gig 2 . ..g k -\llgk+ \ - --gN = gv.k-i^gk+v.N to indicate a 
genotype where g^ = 11 while the background genotype is arbi- 
trary. We will also use the compressed notation 1 IgW (or gener- 
ally gkg (k) )- 

We use the 2-allele content (i.e., the number of 2-alleles) of 
genotypes to define a partial order on the genotype space T (see 
Figure 1, left panel for an illustration). For a particular locus k we 
order the three genotypes sharing the same background genotype 

gl:k-\gk+\:N as follows, 

gl:k-\Hgk + \:N < gv.k- l^gk+ l:N < g\:k- \^gk+ \:N (1) 

We call this the partial genotype order relative to locus k, and it 
defines a strict partial order on T . 

A genotype-phenotype map is a mapping G that assigns to 
each genotype g € T a real-valued genotypic value G(g) (the 
mean trait value for a given genotype). We define monotonicity of 
G in terms of how it transforms the partial genotype order to the 
algebraic order of the genotypic values G(g). Without loss of gen- 
erality we assume that the allele indexes at each locus have been 
chosen such that G(l 1 1 1 - - - 1 1) is the smallest of all homozygote 
genotypic values. We call a genotype-phenotype map G monotone 
or order-preserving with respect to locus k if it preserves the partial 
genotype order relative to locus k, i.e., if, 

G(gl:k-lHgk+l:N) 5 G(^i;jt_il2^ +1:JV ) 

< G(g l:k -i22g k+l ; N ) (2) 
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FIGURE 1 | Examples of partial genotype order and 
genotype-phenotype maps. Left panel: The allele content defines a partial 
order on genotype space. A two-locus example is shown. The plot at the 
top displays the genotype at locus 1 (x-axis) and locus 2 (color) vs. the total 
number of 2-alleles (y-axis) in the two-locus genotype. The resulting partial 
ordering of genotypes is shown below. Right panel: Each lineplot shows 
the 9 genotypic values (y-axis) for a single GP map, coding of genotype are 
the same as in the left panel. GP maps that preserve the partial order of 



genotypes are called monotone. Examples shown are an intra- and 
interlocus additive map (A), a map showing partial dominance at both loci 
(PD), and duplicate dominant (DD) epistasis (see Table 1 in Phillips, 1998). 
GP maps that break the partial order of genotypes are called non-monotone, 
examples shown are pure overdominance at both loci (OD), 
additive-by-additive epistasis (A x A) and dominance-by-dominance epistasis 
(D x D). The rightmost plot shows a GP map that is monotone w.r.t. locus 
1, but non-monotone w.r.t. locus 2. 
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for all genetic backgrounds of locus k. By allowing non-strict 
inequalities we include GP maps showing complete dominance 
and complete magnitude epistasis (Weinreich et al, 2005) in 
the class of order-preserving GP maps. Conversely we call a GP 
map non-monotone or order-breaking with respect to locus k if it 
does not preserve the partial genotype order relative to locus k 
for all backgrounds. Figure 1 (right panel) shows classical dom- 
inance and epistasis patterns, categorized into monotone and 
non-monotone GP maps. 

STATISTICAL DECOMPOSITION OF GENOTYPE-PHENOTYPE MAPS 

Given a genotype-phenotype map G as described above and 
a corresponding vector of genotype frequencies / in a pop- 
ulation, quantitative genetic provides methods for orthogonal 
decomposition of genotypic values and resulting genetic vari- 
ance in the population into additive and non-additive (dom- 
inance and epistasis) components (Lynch and Walsh, 1998). 
We performed such statistical decomposition with the func- 
tion linearGPmapanalysis in the R package noia (http:// 
cran.r-project.org/package=noia; Le Rouzic and Alvarez-Castro, 
2008) version 0.94.1. We assumed an idealized population where 
all genotype frequencies are equal (1/3 N ). In such a hypothetical 
population the NOIA (Alvarez-Castro and Carlborg, 2007) statis- 
tical and functional formulations and the unweighted regression 
model proposed by Cheverud and Routman (1995) are equiv- 
alent. Furthermore, the decomposition of genotypic values is 
equivalent to decomposing G into a sum of additive and non- 
additive GP maps, and the genetic variance in this case is simply 
the variance of the 3 N genotypic values in G. We used the NOIA 
statistical formulation to decompose a GP map G into its addi- 
tive and non-additive components, and computed the ratio of 
additive to total genetic variance Va/Vq as a measure of how 
well the additive component describes the original GP map. In 
case of the illustrative GP maps depicted in Figure 1, this gives 
V A /V G = 1 for the fully additive GP map A, and V A /V G = 0 for 
the pure overdominance (OD) and the pure epistasis (Cheverud 
and Routman, 1996) maps A x A and D x D. 

GENE REGULATORY NETWORK MODELS 

Gene expression in eukaryotes is controlled through gene regu- 
latory networks involving numerous regulatory mechanisms [see 
e.g., Latchman (2005), for details]. Modeling of such gene regula- 
tory networks is well-established, and available modeling frame- 
works range from coarse-grained descriptions of the topology 
of genome-wide networks to very detailed mechanistic models 
describing the dynamics of small networks (De Jong, 2002; Schlitt 
and Brazma, 2007; Karlebach and Shamir, 2008). In line with 
a large number of authors we used ordinary differential equa- 
tions (ODEs) to study a family of generic gene regulatory network 
models containing three diploid genes X\, X2, andX^, organized 
as a regulatory system where the rate of expression of each gene 
can be regulated by the expression level of one or both of the other 
genes. The wiring of the system is described by a 3 x 3 connec- 
tivity matrix A with elements Ay € {—1,0, 1}. The signs of the 
elements of A describe the mode of regulation, Ay = 0 indicates 
that X[ is not a regulator of X^, if Ay = 1 then X/ is an activa- 
tor of Xk, and if Ay = — 1 then Xi is a repressor of X^. Gene 



regulatory systems are often laid out visually as signed directed 
graphs. There is a one-to-one correspondence between a con- 
nectivity matrix and a signed directed graph, two examples are 
illustrated in Figure 4. We used the sigmoid formalism (Mestl 
et al, 1995; Plahte et al, 1998) in the diploid form (Omholt et al, 
2000) where the expression the two alleles of gene k is described 
by the following ODEs, 

Xiu = ukiRki(yuyi,y3) - y*i**i, (3) 
Xk2 = oii2Rjt2(yi, yi,y-i) -YkiXki, 

yk=x k i+x k2 , k= 1,2,3. 

Here ajy is the maximal production rate for allele i of gene X k , 
yid is the decay rate, while Rki is the gene regulation function 
(dose-response function). IfXjfc has no regulators, we assume pro- 
duction is always switched on i.e., Rb= 1. If Xk has a single 
regulator X/, the gene regulation function is given as Rkiiyi) = 
S(yi, Oflti, piki), where S(y, 6, p) = + 6?) if X; is an activa- 

tor and S(y, 6, p) = 1 — y? /(y p + 9?) if it is a repressor. In both 
cases the parameter Gjy gives the amount of regulator needed 
to get 50% of maximal production rate, and p/fa determines the 
steepness of the response. In the case of two regulators X/ and 
Xj we set Rkiiyu Vj) = S(yi, Qm, Piki)S(y P Qjk„ pjki)> corresponding 
to the Boolean AND function. Modeling transcription regulation 
by means of Hill functions and Boolean composition has a long 
tradition in modeling of gene regulation and is widely used. 

With three genes and up to two regulators per gene the number 
of possible connectivity matrices is 6859. We further required that 
the system is connected, and that X3 is downstream to both X\ 
and Xi so either Xi and X2 both regulate X3 directly (A31A32 ^= 
0), or one of them regulates X3 directly and the other one indi- 
rectly (A31A12 ^ 0 or A32A21 5^ 0). This reduces the number of 
distinct connectivity matrices to 3724. Finally, we identified pairs 
of matrices that are symmetric with respect to interchanging Xi 
and X2 and picked just one matrix from each pair. The result- 
ing 1881 connectivity matrices were used for our gene regulatory 
simulations. 

IDENTIFYING FEEDBACK LOOPS AND FEEDFORWARD MOTIFS 

Feedback and feedforward motifs appear recurrently as regu- 
latory building blocks in transcription networks across all liv- 
ing organisms. These network motifs have several characteris- 
tic features (Alon, 2007), negative feedback can for example 
accommodate fast transcriptional responses and homeostasis, 
while positive feedbacks are utilized as biological switches. We 
went through all 1881 gene regulatory models and extracted 
information about their feedback and feedforward loop charac- 
teristics from their connectivity matrices. For each system we 
computed three autoregulatory feedback loop products FL\ = 
An, FL2 = A22, FLt, = A33, three two-gene feedback loop prod- 
ucts: -FI12 = A21A12, FI13 = A31A13, -Fl23=A 2 3A 3 2 and two 
three-gene feedback loop products: FI123 = A32A21A13, FL213 = 
A31A12A23. Non-zero loop products indicate that the system con- 
tains the corresponding feedback loop, and the sign of the loop 
product gives the sign of the feedback loop. We also computed 
the products for two feedforward motifs: FFI32 = A32(A3iAi2), 
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FFL31 = A31XA32A21). Again non-zero products indicate that the 
system contains the corresponding feedforward motif, a positive 
value corresponds to a coherent feedforward while a negative 
value indicates incoherent feedforward. Figure 4 depicts the con- 
nectivity matrix and the signed digraphs of a system with a 
positive feedback loop as well as a system with incoherent feed- 
forward. Spreadsheet SI contains adjacency matrices and loop 
products for all 1881 motifs. 

GENE REGULATORY NETWORK SIMULATIONS 

The simulation were performed with the Python package 
cgptoolbox (http://github.com/jonovik/cgptoolbox), using 
the sigmoidmodel submodule, which contains an implemen- 
tation of the gene regulatory network model (Equation 3) and 
the connectivity matrix A. A similar simulation setup is found 
in Gjuvsland et al. (2011) together with a discussion of gene 
regulation functions and the genotype-parameter map in molec- 
ular terms. We compared two different types of genotype-to- 
parameter maps: 

• Genotype to parameter map without pleiotropy: biallelic geno- 
typic variation for all three loci was introduced through 
the maximal production rates otjt,. For each Monte Carlo 
simulation the allelic parameter values were sampled from 
[7(100, 200). 

• Genotype to parameter map with pleiotropy: allelic parameter 
values were sampled for maximal production rates (sam- 
pled from [7(100, 200)), regulation thresholds Qiki (sampled 
from [7(20, 40)), and regulation steepnesses pm (sampled from 
[7(1, 10)). 

All decay rates yki were set equal to 10. We assembled param- 
eter sets for all 27 diploid genotypes, and for each genotypic 
parameter set the system of Equation 3 was integrated numer- 
ically until convergence to a stable state. The equilibrium value 
of y3 was recorded as phenotype. Datasets where the system 
failed to converge for one or more genotypes were discarded. 
For each of the 1881 motifs we performed 1000 Monte Carlo 
simulations. 

Some Monte Carlo simulations lead to very little phenotypic 
variation, in the sense that the span between the largest and small- 
est of the 27 genotypic values was small. In order to avoid artifacts 
arising from the numeric ODE solver tolerance, these essentially 
flat GP maps were discarded. Further analysis of monotonic- 
ity and variance components were only performed on GP maps 
where the absolute range (maximum genotypic value - mini- 
mum genotypic value) and relative range (absolute range/mean 
genotypic value) were both > 0.01. 

RESULTS 

MEASURING MONOTONICITY OF GP MAPS 

In the following we present two numerical measures for quan- 
tifying monotonicity in a GP map G with N biallelic loci. The 
first quantifies the monotonicity for individual loci by comparing 
negative and positive allele substitution effects before weighting 
the individual loci into an overall measure. The second utilizes 
isotonic regression to quantify the distance between G and the 
closest fully monotone GP map. 



Measure 1: quantifying non-monotonicity by substitution effects 

We first develop a measure of monotonicity based on the effects 
of substituting a single allele at locus k, 

s 1 ^)) = G(g 1:jt _ 1 22 a+1:N ) - Gfcnl2 ft+I;N ), (4) 
s 2 (g (k) ) = Giguk-illgk+ux) - Gduk-illgk+uu), 



while keeping the background genotype g^ = gi-.k+igk+i-.N 
fixed. Monotonicity as defined by Equation 2 is equivalent to 
s'(gW) > 0 for i = 1, 2 across all genetic backgrounds of locus 
k. By taking into account also the magnitude of the substitution 
effects we can quantify the deviation from strict monotonicity. We 
start with the set S k = {s'(g^}} of single allele substitution effects 
for locus k for i = 1,2 and across all genotypic backgrounds gW. 
The total number of elements in S k thus becomes 2 • 3 N ~ 1 , and we 
split the set into two disjoint subsets reflecting their sign; S*j_ = 
{s'(g^) e S k \s , (g^) > 0} and S k i = W(g (k) ) 6 S k \s , (g {k) ) < 0}. 
We compute the sum of positive substitution effects and the sum 
of absolute values of negative substitution effects, 



s'(g (k> ), 



(5) 



and let T k = P k + Nfc denote the overall sum of absolute substi- 
tution effects. We then define the degree to which the GP map G 
is monotone with respect to locus k by, 



m k 



\Pk~N k \ 



(6) 



The absolute value in the numerator ensures that the measure m k 
is invariant with respect to the choice of indexes for the two alle- 
les of locus k. Interchanging the numbering of the alleles leads to 
the mappings s 1 ^)) ^ -j^W), s 2 (g (k) ) \-y s 1 (g^), which 
leaves the value of m k unchanged. By the triangle inequality 
mic < 1. If mfc = 1, then G is monotonic with respect to locus k, 
whereas m k < 1 implies that G is order-breaking w.r.t. locus k. If 
m k = 0, then the positive substitution effects equal the negative 
substitution effects in magnitude and we say that G is completely 
order-breaking w.r.t. locus k. This measure distinguishes well 
between the monotone and non-monotone maps in Figure 1. 
Clearly m.\=mi = \ for the additive map (A) and GP maps 
showing partial dominance and duplicate dominance epistasis. In 
contrast, m\ = m^ = 0 for the maps showing pure OD and pure 
epistasis (A x A and D x D). 

In order to quantify the overall monotonicity of the GP map G 
we introduce the degree of monotonicity (m) which is a weighted 
mean of all m k , where the weights reflect the relative effect size of 
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the loci in terms of Tjfe, 

N 

E m k T k 



E n 

k=l 

As shown in Figure 3A, the degree of monotonicity is accordingly 

I for the monotone maps in Figure 1 while it is 0 for the pure OD 
and pure epistasis maps. This definition of degree of monotonic- 
ity allows us to establish a vocabulary that is analogous to the 
classification of single locus dominance; i.e., a GP map is called 
monotone if m = 1, (partially) non-monotone if m < 1 and purely 
non-monotone if m = 0. 

For example, the degree of monotonicity of the GP map 
published by Cheverud and Routman (1995), with two loci 
underlying 10-week body-weight (in grams) at 10 weeks in a 
mouse F2 cross, may be computed as follows. After renaming the 
two loci (B — >■ 1, A — >2) and indexing alleles to conform to our 
notation, the nine genotypic values (Table 1 in (Cheverud and 
Routman, 1995)) are G(llll) = 31.23, G(1112) = 34.13, 
G(1122) = 33.82, G(1211) = 34.89, G(1212) = 35.90, 
G(1222) = 36.53, G(2211) = 34.12, G(2212) = 37.95, and 
G(2222) = 36.84. From the line plot of this GP map (Figure 2, 
left panel) we find that the GP map is non-monotone with 
respect to both loci. Locus 1 shows marginal OD for the 11 
genotype of locus 2 and locus 2 shows marginal OD for the 

II and 22 genotypes of locus 1. To compute the degree of 
monotonicity, we start with the set of single allele substitution 
effects for locus 1, S 1 = {3.66, -0.77, 1.77, 2.05, 2.71, 0.31}, 
and divide this into sets of negative s! = {—0.77} and pos- 
itive effects S\ = {3.66, 1.77, 2.05, 2.71, 0.31}. The sum Ni 
of elements in S\_ is 10.50 and Pi the sum of absolute values 
of elements in Si is 0.77, which gives T\ = Pi + N\ = 11.27. 
From Equation 6 it follows that mi = 0.86. Similarly, the sets 
of substitution effects for locus 2 are S_ = {—1.11, —0.31} and 
S 2 + = {3.83, 0.63, 1.01, 2.90}. This gives, N 2 = 1.42, P 2 = 8.37, 
T2 = 9.79, and m 2 = 0.71. Inserting values for both loci into 
Equation 7, the degree of monotonicity (m) of this GP map 
is calculated to be 0.79. This value concords well with the 



visual observation (Figure 2, left panel) that it does not deviate 
substantially from a purely monotone map. 

For random GP maps (randomly sampled genotypic values 
as in (Gjuvsland et al., 2011)) there is a strong positive corre- 
lation between the degree of monotonicity and the size of the 
additive component ( Va/Vg) (Figure 3A). A similar relationship 
was observed for three-locus random GP maps (Figure Al A). 
All GP maps in Figure 3A with m < 0.1 have Va/Vg < 0.1. At 
the other end of the spectrum there is much more variation, 
for instance the most extreme completely monotone map (the 
duplicate dominant factors DD) has Va/Vg as low as 0.375. 

Measure 2: quantifying monotonicity by isotonic regression 

This measure quantifies the monotonicity of a particular GP map 
G in terms of the least-squares distance to the closest monotone 
map. We build on the mathematical notation introduced in sec- 
tion "Background on monotonicity of GP maps" where T is the 
genotype space for N biallelic loci and a GP map is a function 
that assigns a real-valued genotypic value G(g) to each genotype 




0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 



v A /v G v A /v G 

FIGURE 3 I Measures of monotonicity vs. additivity of GP maps. 

Scatterplots showing Va/Vq from unweighted regression vs. (A) degree of 
monotonicity (m) and (B) R^om f rom isotonic regression. Black dots 
correspond to the maps shown in Figure 1 together with 
additive-by-dominance epistasis (A x D), a map with two loci showing 
complete dominance (CD) and two classical epistasis types from Table 1 in 
Phillips (1998); duplicate recessive genes (DR) and recessive epistasis (RE). 
Red dots show 1000 random two-locus GP maps, while blue dots show the 
same 1000 GP maps after rearranging genotypic values to introduce 
order-preservation for 1 locus [see Model and Methods in Gjuvsland et al. 
(2011)]. 
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— 22 

— 12 
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11 
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FIGURE 2 I Decomposition of genotype-phenotye map into 
monotone and non-monotone components. Left panel: 

Genotype-phenotype map G for two loci underlying 10-week 
body-weight at 10 weeks in a mouse F2 cross. The GP map shown 
here is equivalent to the one in the original publication [see Figure 



3A in Cheverud and Routman (1995)], but we have changed indexing 
of loci and alleles for consistency with the notation used here. The 
GP map G is decomposed with isotonic regression into a (middle 
panel) monotone component Gm and a (right panel) non-monotone 
component G N . 
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g in T. For any particular GP map G, we identify the monotone 
component of G as the map Gm which minimizes the residual vari- 
ance var(G — Gm), i-e-, Gm is the monotone GP map which is 
closest to G in the least-squares sense. For a given G the mono- 
tone component Gm is unique (Barlow and Brunk, 1972) and 
can be computed numerically by isotonic regression (Leeuw et al., 
2009) of G subject to the partial ordering of genotypes defined 
in Equation 1. Furthermore, the residual Gn is orthogonal to 
Gm in the sense that J^ger GmQOGjvQ?) = 0. This allows the 
orthogonal decomposition, 

G = Gm + Gn, (8) 

of a genotype-phenotype map into a monotone component 
Gm and a non-monotone component Gm such that var(G) = 
var(GM) + var(Gw). The orthogonality property allows us to 
measure monotonicity of G in terms of the coefficient of deter- 
mination i^ ono of the isotonic regression given by the ratio 
^mono = var(GM)/var(G). In the case that G itself is monotone 
for all loci we have J^ono = x » while order-breaking for one or 
more loci will result in _R^ ono < 1. 

The isotonic regression approach can be illustrated in a 
straightforward way on the two-locus GP map provided by 
Cheverud and Routman (1995) (see text above and left panel of 
Figure 2). The partial ordering of genotypes defined by Equation 
1 is illustrated in Figure 1 (left panel). By isotone regression 
(Leeuw et al., 2009) on this partial genotype ordering, the original 
GP map is decomposed into a monotone and a non-monotone 
component (Figure 2, middle and right panels), and the coeffi- 
cient of determination (B^ ono ) is 0.97. 

Our simulation results for random GP maps show that -R^ono 
is positively correlated to the size of the additive component 
(Figure 3B for two-locus GPs maps and Figure A1B for three- 
locus GP maps) and that for a given Va/Vg the lower bound for 
.Rmono is close to a straight line from (0, 0.2) to (1, 1). However, 
due to the search for the closest monotone GP map, -R^ ono will 
not become zero even for purely overdominant or purely epistatic 
maps. As shown in Figure A2, the two monotonicity measures are 
highly correlated. 

An R package for studying monotonicity in GP maps 

We developed an R package gpmap for studying functional prop- 
erties of GP maps. The package takes GP maps in the form of 
vectors of genotypic values as input, and provides functions for 
(i) determining whether the map is order-breaking or order- 
preserving w.r.t. any given locus, (ii) the degree of monotonicity 
m, (iii) ^ mon0 using isotonic regression from the isotone pack- 
age (Leeuw et al, 2009), and (iv) plots of the original and 
decomposed GP maps. Code example 1 (Box 1) below illustrates 
the usage and functionality of the gpmap package. The package is 
available from CRAN http://cran.r-project.org/package=gpmap 
under GPLv3. 

MONOTONICITY IN GP MAPS ARISING FROM GENE REGULATORY 
NETWORKS 

To search for generic relationships between monotonicity and 
regulatory network structure, we used the above measures of 



monotonicity to characterize GP maps emerging from the gene 
regulatory network models (see Models and Methods). Based on 
earlier results (Gjuvsland et al, 2007, 2011; Wang et al., 2013) we 
hypothesized that incoherent feed forward (Figure 4, right panel) 
or positive feedback (Figure 4, left panel) would be necessary in 
order to obtain highly order-breaking GP maps, and we charac- 
terized all 1881 networks in terms of these two properties. Table 1 
shows the number of motifs falling into the resulting four cate- 
gories. We summarized the number of Monte Carlo simulations 
where all genotypic parameter sets gave convergence to a stable 
steady state, and where the resulting GP maps were not essentially 
flat (see Models and Methods for details). Motifs with less than 
100 usable GP maps were discarded from further analysis. For 
the genotype-to-parameter maps without pleiotropy (in the sense 



Box 1 | Code example 1. 

Code example for quantifying and visualizing monotonicity for 
the two-locus GP map published in [14] using the R package 
gpmap . 

> library (gpmap) #load package 

> data(GPmaps) #load dataset 

> gp <- mouseweight #GP map from reference 
[14] 

> 

> ## Tabulate genotypic values 

> cbind (gp$genotype, gp$values) 

> 

> ## Plot the GP map 

> plot(gp) 

> 

> ## Compute degree of monotonicity 

> gp <- degree_of_monotonicity (gp) 

> gp$degree .monotonicity . locus 

> print (gp) 

> 

> ## Quantify monotonicity by isotonic 
regression 

> gp <- decompose_monotone (gp) 

> print (gp) 

> 

> ## Plot decomposed GP map 

> plot (gp, decomposed=TRUE) 
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FIGURE 4 | Connectivity matrices and signed directed graphs. 

Connectivity matrix A and the corresponding signed directed graph for two 
of the 1881 systems in the simulation study. The left panel depicts the 
connectivity matrix and the signed digraph of a system with a positive 
feedback loop between X-\ and X2 while the right panel shows a system 
with incoherent feedforward from to X3. 
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Table 1 | Frequencies (proportion of row total in parenthesis) of incoherent feedforward and positive feedback loops in subsets of the 1881 
studied motifs. 



Dataset Number of motifs 




Motifs containing 




Incoh. feedforward 




No incoh 


feedforward 


Positive feedback 


No positive feedback 


Positive feedback 


No positive feedback 


All motifs 1881 287(0.153) 


48 (0.026) 




1294 (0.688) 


252 (0.134) 


GENOTYPE-TO-PARAMETER MAP WITHOUT PLEIOTROPY 


Discarded motifs 868 152(0.175) 
Analyzed motifs 1013 135(0.133) 


0 

48 (0.047) 




715 (0.824) 
579 (0.571) 


1 (0.001) 
251 (0.248) 


GENOTYPE-TO-PARAMETER MAP WITH PLEIOTROPY 


Discarded motifs 791 124(0.157) 
Analyzed motifs 1090 163(0.149) 


0 

48 (0.044) 




667 (0.84) 
627 (0.575) 


0 

252 (0.231) 



A B 




0 1 2 3 0 1 2 3 

Number of loci for which the GP map is order-breaking 

Fn Incoherent FeedForward & Positive feedback No incoherent FeedForward & No positive Feedback 

P3 Incoherent FeedForward & No positive Feedback c±rj No incoherent FeedForward & Positive Feedback 



FIGURE 5 | Order-breaking in motifs containing a single feedforward 
loop. Summary of order-breaking for all motifs for which at least 100 (out of 
1000) Monte Carlo simulations lead to GP maps with non-negligible 
variation (see Models and Methods section "Gene regulatory network 
simulations," for detailed criteria). Results are shown for 1013 motifs with a 
genotype-to-parameter map without pleiotropy (A) and 1090 motifs with a 
genotype-to-parameter map with pleiotropy (B). Colors indicate classes of 
motifs based on the presence/absence of incoherent feedforward and 
positive feedback loops, see Table 1 for the number of motifs in each class. 
A single boxplot summarizes, for all motifs in the given class, the proportion 
of the GP maps (y-axis) that are order-breaking with respect to a given 
number of loci (x-axis). For example, consider the red box at x = 0 in panel 
(A). This boxplot contains results for motifs with both incoherent 
feedforward and positive feedback and from Table 1 we find that the red 
boxplot summarizes results for 135 motifs. From the y-axis we find that at 
least half (box median at y = 1 ) of these 135 motifs result in only monotone 
GP maps, while for the most extreme (end of whisker) of the 135 motifs 
only 25% of the GP maps are monotone. Similarly, the cyan box is 
compressed into a line at x = 0, y = 1 indicating that all 251 motifs that 
lack both incoherent feedforward and positive feedback result in only 
monotone GP maps. 



that genetic variation at one locus influences only a single param- 
eter, see Model and Methods) 868 motifs were discarded, while 
for the genotype-to-parameter map with pleiotropy (genetic vari- 
ation at one locus influences three parameters) 791 motifs were 
discarded. All (but one) discarded motifs contained at least one 
positive feedback loop (Table 1). A plausible explanation for this 
is that many motifs with positive feedback loops have a stable 
steady state at, or very close to 0 for one or more state variables 
regardless of parameter values, and this leads to essentially flat GP 
maps. 

The introduction of pleiotropy in the genotype to parame- 
ter map has a marked effect on the monotonicity characteristics 
of the associated GP map (Figure 5). When genetic variation 
at a locus X; affects only its maximal production rate the GP 
maps come out as highly monotone (Figure 5A), with a large 
majority being fully monotone or order-breaking for just a single 
locus. When genetic variation at locus X, affects the threshold and 
steepness of the dose-response curve in addition to the maximal 
production rate (pleiotropy in the genotype-to-parameter map), 
the majority of GP maps still show order-breaking either for no 
loci or just one locus (Figure 5B). But a considerable number of 
GP maps are in this case order-breaking for two or three loci. 
Furthermore, by dividing the motifs into the four groups given 
in Table 1 it is evident that the regulatory anatomy of a network 
determines its predisposition for non-monotonicity in its asso- 
ciated GP map. Presence of incoherent feedforward or positive 
feedback loops appears to be prerequisites for the majority of the 
observed non-monotonic GP maps. 

The class of motifs lacking both incoherent feedforward and 
positive feedback contains very few order-breaking GP maps, 
and with no pleiotropy in the genotype-to-parameter map we 
observe only fully order-preserving GP maps for this class (cyan 
in Figure 5A). In the Appendix we generalize this to an arbitrary 
number of nodes and formally prove that without pleiotropy in 
the genotype-to-parameter map, the presence of incoherent feed- 
forward or positive feedback is indeed a necessary condition for 
non-monotone GP maps to arise from networks with monotone 
gene regulation functions. 

The introduction of pleiotropy in the genotype-to-parameter 
map increases the frequency of order-breaking GP maps substan- 
tially (Figure 5B). Motifs lacking both incoherent feedforward 



and positive feedback may in this case lead to GP maps that are 
order-breaking for one or two loci, but never for all three loci. 
Using isotonic regression to quantify the overall monotonicity 
of the GP maps reinforces the finding that incoherent feedfor- 
ward and positive feedback predispose for non-monotonicity 
(Figure 6). Figure 6 also shows that for all classes of motifs 
the majority of GP maps are fully monotone, while the most 
non-monotone GP maps (lowest _R^ ono values) are observed for 
motifs with positive feedback. The differences between classes 
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FIGURE 6 | Empirical distribution functions for R^, ono . Summary of 
^mono values from isotone regression for all motifs for which at least 
100 (out of 1000) Monte Carlo simulations lead to GP maps with 
non-negligible phenotypic variation (see Models and Methods section 
"Gene regulatory network simulations," for detailed criteria). Results are 
shown for 1013 motifs with a genotype-to-parameter map without 



pleiotropy (A) and 1090 motifs with a genotype-to-parameter map with 
pleiotropy (B). Each panel is divided into 4 subplots containing classes 
of motifs based on the presence/absence of incoherent feedforward 
and positive feedback loops, see Table 1 for the number of motifs in 
each class. Each curve shows, for a single motif, the empirical 
distribution function value (y-axis) of ft^ ono for all GP maps (x-axis). 



of motifs are also evident when inspecting the additivity of 
GP maps (Figure A3), but since monotone GP maps can still 
be non-additive, the patterns are much more blurred than for 
monotonicity. 

DISCUSSION 

Fisher's (1918) regression on gene content and the concepts 
derived from this, such as additive effects and dominance devia- 
tion, provide the theoretical basis for most of quantitative genetics 
(Falconer and Mackay, 1996; Lynch and Walsh, 1998). By regress- 
ing on gene content, including the extensions by Cockerham 
(1954), the genotype-phenotype map is decomposed into addi- 
tive, dominant, and epistatic components. The use of gene con- 
tent or the number (0, 1, or 2) of alleles with a particular index 
in a genotype implies the same partial ordering of genotype 
space as defined in Equation 1. Thus, our proposed definition 
of monotonicity of GP maps, and in particular the use of iso- 
tonic regression to quantify monotonicity, may be viewed as a 
relaxation of the linearity assumption underlying current quan- 
titative genetics theory. In this perspective the positive correlation 
between monotonicity and additivity (Figure 3) is expected. 

We have addressed GP maps with 2 and 3 loci as we consid- 
ered an in-depth study of the properties of GP maps with higher 
number of loci to be outside the scope of this study. Some general 
observations can be made, though. Since m is a weighted aver- 
age, the nik of major loci (i.e., for which Xj- is large relative to 
^ T;t) will tend to dominate. For instance, in a case with a single 
major locus showing monotone gene action and several minor 
loci showing order-breaking, the GP map will overall be close 
to monotone (m close to 1). Conversely, order-preservation in a 
number of minor loci would have little influence on m if major 
loci have strongly non-monotone gene action. Isotonic regres- 
sion gives an overall measure of monotonicity of a GP map, but 
provides no locus-specific measures corresponding to m^. Similar 



to the case for m, the gene action of major loci will have high 
influence on the value of R^ono- 

The observation that monotonicity is an important prop- 
erty of GP maps is in principle not new. For a single locus, 
non-monotone gene action appears in the form of over- or 
under-dominance, while complete and partial dominance as well 
as additivity exemplify monotone gene action. Weinreich et al. 
(2005) distinguished between sign epistasis and magnitude epista- 
sis and showed that sign epistasis limits the number of mutational 
trajectories to higher fitness. As sign epistasis reflects a non- 
monotone GP relationship and magnitude epistasis reflects a 
monotone one, this insight concords with our results. A similar 
distinction has been proposed (Wang et al, 2010) for statisti- 
cal interactions where removable interactions are those that can 
be removed by a monotone transformation of the phenotype 
scale, while non-mono tonicity in the GP map leads to essential 
interactions. Wu et al. (2009) developed a method to screen for 
and test the significance of essential interaction in genome-wide 
association studies. Isotonic regression has also recently been 
applied to link genotype and phenotype data (Beerenwinkel et al., 
2011; Luss et al., 2012). Our treatment of monotonicity is more 
general than these earlier works in three major ways. First, we 
deal with monotonicity of the GP map as a whole rather than 
either intra-locus (dominance vs. overdominance) or inter-locus 
(magnitude vs. sign epistasis and removable vs. essential inter- 
actions). Second, where the earlier treatments have focused on 
classifying the type of gene action, we make use of quantitative 
measures of monotonicity. Third, our approach combining the 
concept of monotonicity with cGP models opens a direct link 
between genetics and the theory of dynamical systems in the wide 
sense. 

Monotonicity is a property of the GP map separate from 
the allele frequencies, making it a physiological (Cheverud and 
Routman, 1995) or functional (Hansen and Wagner, 2001) 
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descriptor rather than a statistical one. The distinction between 
physiological and statistical epistasis has lead to much debate 
(Phillips, 2008). Zeng et al. (2005) argued the distinction was 
unnecessary and potentially misleading. Although their argu- 
ments around orthogonality and variance components are valid, 
our results demonstrate very clearly that describing the properties 
of the GP map without reference to any particular study popula- 
tion is essential if we want to connect quantitative genetics with 
regulatory biology. 

It is clear from our results that positive feedback and incoher- 
ent feedforward promote non-monotonicity. The clear-cut dif- 
ferences in monotonicity between different classes of regulatory 
networks, combined with the strong correlation between mono- 
tonicity and additivity of GP maps, appear therefore to explain the 
findings that regulatory systems with positive feedback give con- 
siderably more statistical epistasis than those without (Gjuvsland 
et al, 2007; Wang et al, 2013). Even though both incoherent feed- 
forward and positive feedback predispose for non-monotone GP 
maps, the underlying mechanisms are different for the two reg- 
ulatory motifs. In the case of incoherent feedforward the sum of 
direct and indirect effects may result in a non-monotone dose- 
response relationship (Kaplan et al., 2008). That positive feedback 
loops can give non-monotonicity is intuitively less clear, but in 
the Appendix we show both results analytically. Positive feedback 
predisposes for multiple steady states, and order-breaking might 
also emerge from different genotypes corresponding to different 
states. It should be noted, however that positive feedback is only 
a necessary condition for multistationarity (Plahte et al., 1995), 
and a positive loop in the connectivity matrix A of a system is 
not necessarily active at any point during the time course of the 
system. 



Without any restrictions on the connectivity of a three- 
gene system there are 3 9 = 19, 683 possible distinct networks. 
The main restriction we imposed (see Models and Methods 
for details) was a maximum of two regulators per gene, which 
allowed us to use Boolean gene regulation functions already 
established in the sigmoid formalism (Plahte et al., 1998). Other 
model formalisms allowing an arbitrary number of regulators 
are also available (Wagner, 1994, 1996; Siegal and Bergman, 
2002) and could be extended to diploid forms and used in later 
studies. 

Although this study has focused on gene regulatory net- 
works, the concept of monotone gene action applies to the 
propagation of genetic variation across the whole physio- 
logical hierarchy. One may therefore systematically use the 
concepts and methods presented here to study the order- 
preserving and order-breaking properties of genotype-phenotype 
mappings that are associated with any regulatory structure 
amenable for mathematical modeling. Through this it will be 
possible to make a wide-ranging survey of which regulatory 
anatomies promote monotonicity and which promote non- 
monotonicity. We foresee that this classification may become 
instrumental for predicting how phenotypic effects of genetic 
variation propagate across generations in sexually reproducing 
populations. 

SUPPLEMENTARY MATERIAL 

The Supplementary Material for this article can be found online 
at: http://www.frontiersin.org/journal/10.3389/fgene.2013.00216/ 
abstract 

Spreadsheet SI | Excel spreadsheet with connectivity matrices and loop 
products for all 1881 gene regulatory networks. 
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APPENDIX 

In this appendix we complement the simulation studies in the 
main text with some analytic results for GP maps emerging from 
ODE models of gene regulatory networks. We study a general- 
ization of the gene network model in Equation (3) with an arbi- 
trary number of loci and monotone gene regulation functions, 
but restrict the analysis to genotype-parameter maps without 
pleiotropy. In particular, we show that (i) if there are no posi- 
tive feedback loops and no incoherent feedforward loops in the 
network, the resulting GP maps are always monotone, (ii) a posi- 
tive feedback loop or an incoherent feedforward loop may lead to 
non-monotone GP maps. The results hold for phenotypes given 
as the stable concentration of the product of one of the genes, 
and under certain restrictions also for phenotypes given as a func- 
tion of one or several stable gene product concentrations that is 
monotonic with respect to each of its arguments. 

GENE NETWORK MODEL 

We consider a dynamic system consisting of n mutually inter- 
acting diploid lociXj, j e N = {1, . . . , n), regulating each other's 
expression. The time dependent output of X; is denoted Zj, and we 
define z = [zi, Zi, . . . , Z«]. It goes without saying that Zj in gen- 
eral depends on the genotypes of all the genes even though we will 
not always state this explicitly. 

For a given genotype g = Kg® = «/&jg® , where 
gj € {11, 12, 22} denotes the genotype and a;, bj G 1,2 denote 
the indexes of the two alleles of locus Xj, the equations of motion 
for X, are 





= <V 


1 a U \ 






b 

= a-' 
) 


'^(Z) 


b i 2 




= 4 







where zj and zj are the time-dependent outputs of the two 
homologous copies of Xj. The two allele rate functions rj (z) and 
r?(z) have range [0, 1] so that aj and aj represent the max- 
imum production rates of the two alleles. We assume that all 
dose-response functions in Equation (Al) are differentiable and 
monotonic with respect to each of its arguments, and that for each 
k, the signs of drj/dxk and drj/dxk in the stable point x are 
equal. This model generalizes Eq. (3) to an arbitrary number of 
loci and a broader class of gene regulation functions. 

In the following we are only concerned with the steady states 
of Equation (Al), and assume for simplicity that they have just 
a single stable equilibrium point. Solving the equilibrium condi- 
tions of Equation (Al) with respect to zj and zj and adding gives 

f j ( X ) = [ L?r't i (x) + v b j i rf(x)-Xj = 0, j e N, (A2) 
where x = [x\, . . . , x n ] is the stable point, [ij' = aj 1 /yj and 

b- b- b 

\h-' = a-' /y-' . Since our definition of monotonicity of GP maps 
does not depend on the numbering of alleles, we will without loss 
of generality assume (i 1 < for all j. 



The network architecture can be read out from the structure 
of the system's Jacobian matrix in the stable state x. We define the 
elements of the Jacobian / for the set of functions fi defined in 
Equation (A2) by 

Ijk = Jjk(g) = ^— , j,ke N. (A3) 

OXk 

To the Jacobian / it is customary to assign a signed directed graph 
Q in which each locus X% is represented by a node Xjt, and in 
which there is an arc from Xj to X^ if and only if Ju / 0, its sign 
given by the sign of Ju. A chain from Xj to Xt is a set of arcs in Q 
leading from Xj to X,t in which all intermediate nodes are visited 
only once. The sign of a chain is equal to the product of the signs 
of the Jij corresponding to the arcs in the chain. If there is a chain 
from X, to Xj and also a chain from Xj to X, through a disjoint set 
of nodes, the two chains constitute a proper feedback loop (FBL). 
To each FBL is associated a loop product L which is the product 
of the Jacobian elements corresponding to all the arcs in the loop. 
The sign of the loop is given by the sign of L. Two chains from Xj 
to X,, i ^ j, with only the endpoint nodes in common, constitute 
a feedforward loop (FFL). If the two chains have opposite signs, 
the FFL is incoherent (IFFL), otherwise it is coherent (CFFL). 

The system's phenotype could be any scalar quantity defined 
by its equilibrium value x. In the following we assume the 
genotype-phenotype map G(g) = x q (g), q e N, for a given and 
fixed q, and investigate the monotonicity properties of G(gkg^) 
with respect to genetic variation in any locus X^ for different back- 
grounds gw . i n the following sections we analyse the causes of 
order-breaking in G in the restricted case in which there is only 
genetic variation in |Xjf and \x? k , not in the shape of the dose- 
response functions rj and rf, implying rl(x) = ri(x) = 
This is what we mean by a genotype-to-parameter map without 
pleoitropy. 

In the next sections we prove the following result: 

Proposition 1. Assume all rate functions in Equation (Al) are 
monotonic and that G is the mapping from g to x q for some fixed q 
so thatx q (g) is the phenotype. If there is no feedback loop (FBL) and 
no feedforward loop (FFL) anywhere in the network corresponding 
to the system Equation (Al ), then necessarily = 1 for all k. If the 
system contains either a single FFL or a single FBL, then G may be 
non-monotone for some xjt if the FFL is positive or the FBL is inco- 
herent, but if the FBL is negative or the FFL is coherent, no order 
breaking can occur for any x\^. 

At the end of this note we show that under some reasonable 
conditions this result is also valid for more general phenotypes 
depending on more than one x q . 

NETWORKS WITHOUT LOOPS 

We first consider networks containing no feedforward loop and 
no feedback loop. In these networks there is at most one chain 
from one node to another, and of course no autoregulatory loops. 
If there is a chain from Xj to X^, there is no chain from X^ to 
Xj. Any node is either unregulated (constitutively expressed) or 
regulated by one or several other nodes. 



www.frontiersin.org 



November 2013 | Volume 4 | Article 216 | 11 



Gjuvsland et al. 



Monotonicity is a key feature of genotype-phenotype maps 



We first prove a useful lemma. 

Lemma 1. If < x/(12gW) < x/(22g©) for any j and I 

and there is an arc X/ — »■ X m with positive sign and no other chain 
from Xi -> X m , then also x m (llg©) < x m (l2g ( ^) < x m (22g^). 
If the sign of the arc is negative, then x m (llg^) > x m (12gW) > 
x m (22g^)- 

Proof. Suppressing the explicit dependence on other genes that 
are not affected by genetic variation in Xj, we have 



x m (llg^) = 2yx, 1 m r m (x l {llg^)), 
x m (12g©) = {\x} m + Li^)r m (x ; (12^'>)), 
X„(22g®) = 2\i 2 m r m (x l (22g<J ) )). 



(A4) 



Now, r m is monotonic by assumption. If it is monotonically 
increasing, 



x m {\2g^) > (m4 + |i£)r m (xj(ll*®)) > XmiUg®), 



x m (22g (i) ) > 2^ m r m (xi(W n )) > XmW) 



Oh 



(A5) 



from which the assertion follows. If r m is monotonically decreas- 
ing, we find the same relations with the inequality signs 
reversed. □ 

If there is no chain from Xj to X q , genetic variation in Xj will 
not be reflected in G, i.e. G(llg®) = G(\2g^) = G(22g©), and 
by definition does not give order-breaking. Then assume Xj is 
upstream relative to X q and that the chain from Xj to X q is posi- 
tive. We first let Xj be an unregulated node with no predecessor. 
Then 



x j {ng®)=2 Vu ), 
Xj(12g®) = u.j + u|, 
^(22^) = 2jjl?, 



(A6) 



because r 



1. From this it follows that Xj{l\g^) 



Xj(UgU) < Xj{22g®), 

Repeated use of Lemma 1 leads eventually to x q (llg^) < 
x q (\2g^) < x q (22gW), irrespective of the genotypic background 
of Xj. If the chain from Xj to X q is negative, the argument goes 
in the same way, but then x q (llg^) > x q {\2g^') > x q (22g^). 
The above argument can be carried out in the same way if Xj 
is not top-stream. It follows that in a network without FFBs 
and FFLs and where genetic variation is restricted to (Jijf and 
(jijr, the genotype-phenotype map G(g) = x q (g) cannot be order- 
breaking. 

NETWORKS WITH A FEEDBACK LOOP 

In this section we investigate the effects of feedback loops on 
the degree of monotonicity. Assuming monotonic dose-response 
functions and non-pleiotropic genetic variation, we show that a 
positive feedback loop may lead to order breaking, while negative 
feedback loops never do. We consider a network in which there is 



no FFL and a single FBL with X q as one of its members and Xk is 
upstream of the loop. 

Lemma 2. Consider a network with n nodes for which all dose- 
response functions are monotonic and there is only genetic variation 
in [il and |xjr. Asssume there is a chain from Xk to X\, that X\, 
but not Xk, is member of a FBL with m nodes, and that there is no 
other FBL and no FFL in the system. If X q is in the loop, let the 



loop be X\ —>■ X2 



Xn 



X m -> Xi. If the FBL is 



positive, there may be order-breaking in X q due to genetic variation 
in Xk, but no order-breaking can occur if the loop is negative. IfX q 
is downstream of the loop, the same result applies. 

Proof. With a single FBL and no FFL there is at most one directed 
path from any node X; to any other node Xj, and if there is a path 
from Xi to Xj, there is no return path from Xj to X; if either X, 
or Xj is not part of the FBL. We first consider the dependence 
of x\ on Xk. The direct regulators of node Xi are X m and X;, the 
latter being the last but one node in the chain from Xt to Xi . In 
Plahte et al. (2013) we introduced the propagation functions Xj = 
Pjk( x k) which express the effect on Xj of genetic variation in X^. 
An important property of pjk is that it can be derived from all 
the equilibrium conditions Equation (A2) except the equation for 
/it. This implies that the effects on Xj of genotypic variation in 
Xjt are only expressed in terms of the variations in xi, while the 
parameters expressing the genotype of Xk do not enter into the 
function pjk . 

We then have x\ = pik(xk) andx,„ = p m i(xi). To make it easier 
to use the results in Plahte et al. (2013) we rewrite the equilibrium 
condition Equation (A2) as 



Rj(x) - YjXj = 0, 



(A7) 



where y; > 0. In the following, the Jacobian refers to this set 
of equations, which has the same root and the same functional 
dependencies between the variables as the original set. The signs 

of the partial derivatives of -R ; - are the same as for rj and r^. The 
equilibrium condition for Xi is then 



Y1X1 = Ri(pik(x k ), p m i(xi))). 



(A8) 



This equation defines x\ as a function of Xk in an open domain 
around the equilibrium point and with a derivative that can be 
computed by implicit differentiation, i.e. 



dxi dR\ dRi dx\ 

Yl-5— = "5 Ilk + qml-7— ■ 

axk ox\ 3x m axk 



(A9) 



where = p'- is the derivative of pij for all i, j. 

From Lemma 1 it follows that there is no order breaking in X;, 
in other words, qik has a fixed sign. Consider then q m \. There is 
just a single chain from Xi to X m , and Equation (13) in Plahte 
etal. (2013) gives 



,i(xi) = (-1)" 



DwCu 



(A10) 



Here U is the set of nodes in this chain, Cjj is its chain prod- 
uct, i.e. the product of the Jacobian elements corresponding to 
the arcs in the chain, V = N \ U, ' is the subdeterminant of / 
with row 1 and column 1 deleted, and Dyv is the subdeterminant 
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of / composed of the rows and columns V. Because there is 
no feedback loop among the nodes represented in D^ n ' and 
Dyv> only the diagonal degradation terms contribute to these 
two determinants. Hence D* 11 * = (— l) n_ 1 Ylijt 1 Yi- Similarly, 

Dw = (-l)" _m n, e y Yi. giving 1ml = YiQ// r u, where = 
rii'euY;- Finally, we note that P= (dRi/dx m )Cu is the loop 
product of the loop. 

Solving Equation (A9) with respect to dx\/dxk and using all 
these expressions lead to 



Yi 



cbci 
dx k 



T v dxi 



P dxi 



qik- 



(All) 



The sign of dxi/dxi is independent of the genotype of Xk and the 
sign of qi^ is fixed. Genotypic variation in Xk may change the 
magnitude of P, but its sign is fixed because all Jacobi elements 
have fixed sign independent of the system parameters. Thus, 
genotypic variation in X k does not alter the sign of dx\ /cbcfc if 
the loop is negative (P < 0), while for a positive loop the sign 
of Fu — P may switch. In the latter case, an increase in x k due to 
genetic variation in Xk may increase x\ in some cases and decrease 
it in others, leading to order breaking. As there is only a single 
chain from X\ to X q , no order breaking in X\ implies no order 
breaking in X q , while order breaking in X\ may propagate to X q . 
The same result follows if X„ is downstream a node in the loop 



because order breaking in this node may propagate to X q 



□ 



FEEDFORWARD LOOPS (FFLS) 

A feedforward loop (FFL) is a motif in the network in which there 
are two different chains C\ and C2 from one particular node to 
another particular node. To each chain C; is associated a chain 
product P[ defined as the product of the Jacobian elements corre- 
sponding to the arcs in Q. If Pi and P2 have equal signs, the FFL 
is coherent, otherwise it is incoherent. 

In a network with a single feedforward loop and no feedback 
loops we now investigate the effect on G(g) = x q (xk(g)) of genetic 
variation in Xk for varying background . Our starting point 
is again Equation (A7). We first let Xk and X q be the initial and 
terminal nodes in the FFL. The two chains Ci and C2 leading 
from Xk to X q comprise pi and P2 nodes including Xk and X q , 
respectively. Let the set of nodes in C\ and C2 be X^ and Xjj 2 > 
respectively, where U\ and U2 are the corresponding subsets of 
N, and let V\ and V2 be their complements. 

Roughly speaking, the derivative of the propagation function 
p q k(xk) can be expressed as a sum of terms, each term correspond- 
ing to one of the chains leading from Xk to X q (Plahte et al., 2013). 
To the chain C, is assigned the chain weight w; given by 



(-l)P-- 



1, 2, 



(A12) 



where Dy^ is the Jacobian subdeterminant for the nodes not 
included in Q, and D^) i s the Jacobian subdeterminant for all 
nodes except Xk. Because there are two chains from Xk to X q , the 
derivative of p q k is a sum of two terms: 



dpqk 
dx k 



W1P1 + W2P2, 



(A13) 



where Pi and P2 are the two chain products, and w\ and W2 their 
weights (Plahte et al., 2013). When there is no feedback loop in 
the system, only the diagonal elements in / stemming from the 
term — y,x,- in Equation (A7) contribute to the determinants Dy i v i 
andD^: 



Dv.v, = (-D"- p 'nY;, 

j€V, 

D (u ) = ( _ 1) „- 1T -[ Y .. 



(A14) 



Altogether this gives 



dxq = dpqk = % p n p 
dx k dx k Ti T 2 



(A15) 



where Ti and T2 are the products of the y ; in the two chains, 
respectively. The chain products Pi and P2 depend on the geno- 
type gk of Xk as well as on the genotypic background g™, but 
their signs Si and S2 are invariant under genotypic variation. It 
is easy to see that a negative autoregulatory loop, which is a com- 
mon feature in gene regulatory networks, would not invalidate 
the conclusion, but a positive autoregulatory loop might. 

If the FFL is incoherent, Pi and P2 have opposite signs, imply- 
ing that the sign of dx q /dxk may vary. If the FFL is coherent, 
however, no order-breaking can occur. 

If Xk is upstream relative to the initial node Xi n ; t of the FFL, 
it follows from the above section on networks without loops that 
there will be no order-breaking in X m i t , and the above argument 
is still valid. 

MORE GENERAL PHENOTYPES 

In real life, relevant phenotypes are not direct gene products, 
but rather functions of the concentrations of one or several gene 
products. Let the phenotype G{g) be a function of x\j{g), G = 
h(xjj(g)), where U is a subset of N, and assume that for any 
u e U, dh/dx u has fixed sign for all genotypes. To analyse this 
case we extend the original system Equation (A2) to 



h b 



Vi'ri'(x(g)) + \l { ' r ^ {x{g)) — x t {g) = 0, i=l,...,n, 



(A16) 



Hxuig)) ~x„ 



+ 1 



0. 



and apply our above results to this system, in which G(g) = 
x n + i(g), i-e- q = n + 1. If there are two nodes among Xjj which 
have a common predecessor Xk, then there will exist two chains 
from Xk to X n + \ . These two chains constitute a feedforward loop 
with X n+ 1 as final node. If this FFL is incoherent, order breaking 
due to genetic variation in Xk may occur even if there is no order 
breaking in the original system comprising the nodes X\ , . . . ,X n . 
If the FFL is coherent, order breaking only occurs if it occurs in 
the original system. 
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FIGURE Al | Measures of monotonicity vs. additivity of GP maps with 
three loci. Scatterplots showing Va/Vq from unweighted regression vs. (A) 
degree of monotonicity (m) and (B) fi^ ono . Red dots show 1000 random 
three-locus GP maps, blue dots show the same 1000 GP maps after sorting 
to introduce order-preservation for 1 locus while green dots show the same 
1000 GP maps after sorting to introduce order-preservation for 2 loci [see 
Model and Methods in Gjuvsland et al. (201 1 )]. 
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FIGURE A2 | Comparing measures of monotonicity GP maps. 

Scatterplots showing degree of monotonicity (m) vs. R^ ono . Black dots 
correspond to the maps shown in Figure 1. Red dots show 1000 random 
two-locus GP maps, while blue dots show the same 1000 GP maps after 
sorting to introduce order-preservation for 1 locus [see Model and Methods 

in Gjuvsland etal. (2011)]. 
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FIGURE A3 | Empirical distribution functions for additivity of GP maps. 

Summary of V A /V G from unweighted regression for all motifs for which at 
least 100 (out of 1000) Monte Carlo simulations lead to GP maps with 
non-negligible phenotypic variation (see Models and Methods section "Gene 
regulatory network simulations," for detailed criteria). Results are shown for 
1013 motifs with a genotype-to-parameter map without pleiotropy (A) 



and1090 motifs with a genotype-to-parameter map with pleiotropy (B). Each 
panel is divided into 4 subplots containing classes of motifs based on the 
presence/absence of incoherent feedforward and positive feedback loops, 
see Table 1 for the number of motifs in each class. Each curve shows, for a 
single motif, the empirical distribution function value (y-axis) of Va/Vq from 
unweighted regression for all GP maps (x-axis). 
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